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An ensemble of uncoupled limit-cycle oscillators receiving common Poisson impulses shows a range 
of non-trivial behavior, from synchronization, desynchronization, to clustering. The group behavior 
that arises in the ensemble can be predicted from the phase response of a single oscillator to a given 
impulsive perturbation. We present a theory based on phase reduction of a jump stochastic process 
describing a Poisson-driven limit-cycle oscillator, and verify the results through numerical simula- 
tions and electric circuit experiments. We also give a geometrical interpretation of the synchronizing 
mechanism, a perturbative expansion to the stationary phase distribution, and the diffusion limit 
! >— ■ ; ■ of our jump stochastic model. 
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I. INTRODUCTION 

The improvement in response reproducibility of a system receiving identical fluctuating drive has received much 
attention recently 0, H |j, 0, [|, B 0, H, IE Eo, HI E3, EE 0, EH- With a constant input signal, many systems 
composed of identical elements in an ensemble show unreliable response due to noise after starting from similar initial 
conditions, while a fluctuating drive greatly improves response reproducibility. This general phenomenon is observed 
in various guises throughout nature. Uncoupled lasers driven by a common master laser with a fluctuating output 
1 ■ produce highly correlated output intensities |l[ . In the rat neocortical neuron [2] , action potential (spike) generation 
times over many trials coincide when a fluctuating current is delivered to the soma. Synchronized firings in vivo in 
the cat spinal motoneurons [|[ , and in neurons of the olfactory bulbs in mice [H have also been reported. Population 
density correlation among spatially separated species, known as the Moran effect is a general phenomenon seen in 
many different organisms in the troposphere. The phenomenon has also been found in chaotic oscillators [E0, where 
00 ' synchronization of the generalized phase has been observed, and is known as noise-induced chaotic synchronization. 
I/"") , Common forcing may also lead to a decrease in response reproducibility [E E2, EH, EE EE EE E3 ■ 
C ' Let us first illustrate with examples the phase coherence phenomenon induced by fluctuating drive. Figures [T] 
and [11] show the orbit of unperturbed limit-cycle oscillators obtained numerically and experimentally, and Figs. [2] 
£NJ ■ and [M] show the effect of common impulses on an ensemble of such oscillators. Though the oscillators do not interact 
with each other, they exhibit phase synchronization and desynchronization [43| depending on the impulse intensity 
00 , and oscillator characteristics. Knowledge of the response of the oscillator to an impulsive perturbation is sufficient to 

' understand the observed coherence phenomena, as we will clarify in this paper. 
J> , Several previous studies have demonstrated this phenomenon for various driving signals and oscillators [IE EH, EE 
■ EE, EE EH by using the phase reduction method for limit cycles [IE [HI- In particular, the case where the fluctuating 
signal is a sequence of random impulses has been investigated in Pikovsky, Rosenblum and Kurths (Q, Sec. 15 and 
references therein) and by ourselves [IE]. Using random phase maps, it is argued that synchronization always occurs 
for general limit cycles when sufficiently weak additive random impulses are given, and desynchronization can also 
occur when the impulse intensity is finite. 

In this article, we generalize our previous argument through a refined formulation in terms of Poisson-driven Markov 
processes [22], also known as jump processes [IE], which enables us to systematically perform phase reduction and 
linear stability analysis of our impulse-driven oscillators for general multiplicative coupling. The oscillator response to 
the impulses is expressed as a function called the phase response curve (PRC) 0, HH , which is a very basic quantity of 
limit-cycle oscillators measured in every experiment, and the coherent behavior that arises from the common impulse 
can be deduced once the PRC is obtained. We present a criteria for predicting when each of the above coherence 
phenomenon can be expected from the application of common impulses, and test the predictions quantitatively using 
numerical simulation and an electric circuit experiment. We measure the PRC for a typical limit-cycle oscillator 
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described by a set of ordinary differential equations, and for an electrical limit-cycle oscillator. We then compare 
the rate of synchronization or desynchronization as predicted by the Lyapunov exponent obtained from the PRC in 
both numerical simulation and experiment. We also give a simple geometrical explanation of the synchronization as a 
consequence of the stability of the limit cycle, a perturbative expansion to the phase distribution of a single oscillator, 
and a derivation of the diffusion limit of the jump process describing our impulse-driven limit-cycle oscillators in the 
Appendix. 

II. THEORY 

In this section, we present linear stability analysis of the impulse-driven oscillator through phase reduction of 
the jump stochastic differential equation describing the model. Similar analyses have been performed based on 
random phase map description for a simple sinusoidal phase map in Ref. Sec. 15, and for more general phase 
maps derived from the phase reduction of limit cycles in Ref. [l3[, which gave formulas relating the Lyapunov 
exponent to the functional form of the phase maps and predicted synchronization and desynchronization of uncoupled 
oscillators driven by common random additive impulses. Our reformulation based on the stochastic jump process 
given in this section provides a simple and systematic treatment of the impulse-driven limit cycle oscillator for general 
multiplicative coupling, which quantitatively relates the Lyapunov exponents, the PRCs, and the phase-space structure 
(isochrons [20] ) in a mathematically transparent way, starting from a general dynamical equation describing impulse- 
driven limit cycles. It also incorporates the effect of different stochastic interpretations for the random impulses. 
Using our new formulation, we derive the classical results in a more general way and argue the possibility of clustered 
states for multiplicative impulses. 

A. Model 

We consider an TV-oscillator ensemble receiving a common sequence of random Poisson impulses. The equation for 
the a-th oscillator is 

N(t) 

X^ a \t) = F(JfW) + <r(X {a \c n )h(t-t n ), (1) 

n=l 

where a = 1, • • • , N, X^(t) £ R M is the oscillator state at time t, F(X ( ")) : R M —> R M the dynamics of a single 
oscillator, N(t) the number of received impulses up to time t, t n the arrival time of the n-th impulse, c n £ R the 
intensity and direction, or mark [12, HH, of the n-th impulse, a(X^ a \ c) : R M x R K — ► R M is the coupling function 
describing the effect of an impulse to X^ a >(t), and h(t — t n ) is the unit impulse (J^° h(t — t n )dt = 1) whose waveform 
is localized at the event time t n . We denote the rate of the Poisson impulses as A, namely, the mean interval between 
the impulses is 1/A. 

In the absence of impulses, we assume that each oscillator obeys the same dynamics, with a single stable limit cycle 
solution Xo(t) of period T in phase space. In the following, we omit the oscillator index a as the ensemble is composed 
of identical uncoupled elements, and our discussion involves only the linear stability of an individual oscillator. 

B. Jump stochastic differential equation 

We pose the problem as a Poisson-driven Markov process, or stochastic jump process. We regard the impulse as an 
event of zero-temporal width, so that the temporal correlation of the impulses vanishes. The resulting discontinuous 
oscillator dynamics can be described by a stochastic- inte gro differential equation, or jump stochastic differential 
equation (jump SDE), for a marked Poisson point process |22l . |23| . Some properties of the jump SDE are given in 
Appendix A. 

There are two ways we may interpret the impulse term in the original ordinary differential equation (fT]), which we 
call Ito and Stratonovich pictures for convenience. The Ito-picture assumes the impulsive term in the original Eq. (fT|) 
as being point-like impulses, which gives rise to discontinuous system dynamics with jumps. The Stratonovich-picture 
assumes the impulsive term in the original Eq. (TT]) as a limit of short but continuous waveforms of non-zero width, 
thus requiring that we find the limit of the system response as the impulses are shrunk to width. 

Both pictures lead to the same jump SDE for the phase variable: 



dX{t) = F(X)dt + / g(X,c)M(dt,dc), 

J c 



(2) 
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which is always interpreted in the Ito sense. Here, M(dt,dc) represents a Poisson random measure, which gives the 
number of incident points during [t, t + dt] having the mark [c, c + dc], whose expectation is 

(M (dt, dc)) = Xdtp(c)dc, (3) 

where A is the rate of the Poisson process and p(c) is the probability density function (PDF) of the marks c, and the 
integral is over the mark space (Ref. [z^. |23|. Appendix A). The oscillator response g(X,c) to a given impulse with 
mark c is different between the two pictures when the effect of the impulse <x(X, c) is multiplicative. When the full 
oscillator dynamics is reduced to phase dynamics, the phase response is slightly altered, as we see later. 



1. Ito picture 

We interpret the common impulse h(t — t n ) in the original Eq. ((T|) to be zero-width from the outset. When an 
impulse c is received at time t, the state of the oscillator jumps discontinuously from X to X + er(X,c). This 
interpretation gives a jump SDE of the form 

dX(t) = F(X)dt + J cr{X,c)M(dt,dc). (4) 

Thus, g(X,c) = a(X,c). 



2. Stratonovich picture 

We interpret the impulse h(t — t n ) in the original Eq. (TTJ) as having a continuous waveform of finite width and then 
shrink the impulse width to 0. The resulting continuous process can be approximated by a discontinuous jump SDE 
with a jump amplitude as determined by the canonical form of the Wong-Zakai theorem for jump processes as shown 
by Marcus [24|, ■ Defining the differential operator 



M 

1=1 



d 
dXi 



(5) 



the jump SDE is given by 



dX(t) = F(X)dt 



e D X-X 



M(dt,dc). 



(6) 



Thus, g(X, c) = e D X — X in this case. 

Note that for additive impulses, <t(X, c) = <r(c), Eq. (TJ|) and Eq. © take the same form because e D X — X = <x(c). 
Therefore, the difference in stochastic interpretation is not important in this case. For the linear multiplicative 
case, <Jk(X,c) = c^Xu- ft is easy to check that gk(X,c) = (e Ck — f)A"fc. This expression for the instantaneous 
jump magnitude approximates the effect of an impulsive perturbation due to a short but finite-width impulse whose 
intensity in each direction is cu- For more general multiplicative coupling, explicit expressions for g(X, c) are difficult 
to obtain. 



C. Phase reduction 

To facilitate theoretical analysis, we apply phase reduction (2^, Hl| to Eq. @, assuming that the average time- 
interval between jump events is long compared to the relaxation time of the perturbation to the limit cycle orbit. 
An unperturbed oscillator executes periodic motion along its limit cycle, so its state can be described by one phase 
variable <fi(t) = (f>(Xo(t)) 6 [0, 1 ) which constantly increases with frequency ui = l/T, instead of the original M 
variables. 

Because the limit cycle is assumed to be globally stable, the orbit of any initial point P off of the limit cycle will 
asymptotically approach the limit cycle. Thus, we can extend the definition of the phase 4> to the whole phase space 
except at phase singular sets by identifying the set of points that asymptotically converge to the same orbit on the 
limit cycle with the same phase, called the isochron [20l. [2f|. In practice, if the orbit approaches a point on the limit 
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cycle with phase <p to within some arbitrarily small distance after time r. we define the asymptotic phase of the initial 
point P to be (</> — t/T) mod 1, where T is the period of the oscillator. 

Now we perform phase reduction, which is mathematically a change of variables describing the system from X 
to <p = cf)(X), and is also an approximation of the function of X by the corresponding function of Xo((j)). In this 
transformation, every value of X in the neighborhood of the limit-cycle attractor maps to a value of <p except at phase 
singular sets. Using the stochastic chain rule for the jump process (Refs. [HIHH, Appendix A), we obtain 

d<j){t) =u;dt+ I [cf) (X + g(X, c)) - 4>{X)} M(dt, dc), (7) 

which is not yet a closed equation for p. The map X — > X + g{X, c) describes the effect of an impulse received when 
an oscillator is at X. Since we assume that the average interval between impulses is longer than the relaxation time 
of the perturbation to the limit cycle, we can evaluate the function of X using values of Xq on the limit cycle, for 
which the mapping <p — > Xq is well defined. Replacing the X with Xo(<f>), we obtain 



d<t>{t) = Ludt + J [4>{X Q {4>) + g{X {4>),c)) - <t,]M{dt,dc) 

= ujdt+ G(<t>,c)M(dt,dc), (8) 



which is now a closed equation for <p. Here we introduced a function 

G(0, c) = 4> (X Q (cf>) + g(X (<f>), c)) - 0, (9) 

which is the phase response curve (PRC) representing the change in asymptotic phase relative to an unperturbed 
oscillator caused by an impulse with mark c received at phase <p. Since we consider a continuous dynamical system, 
G{<p, c) is continuous and periodic in <p. If the jump amplitude g(Xo(p>), c) is small, G(0, c) can be approximated as 

G(0,c)-Z(0).g(X o (0),c), (10) 

where 

Z(0) = grad x 0(X)| x=XoW (11) 

is the well-known phase sensitivity function that represents linear sensitivity of the phase to infinitesimal perturba- 
tions [IlHHl- In general, this linear relationship holds only for very weak impulses (see e.g. [HI). The PRC easily 
becomes a largely fluctuating function, which can even become jagged, e.g. near the bifurcation point. 

Note that for our phase-reduction analysis of Poisson-driven limit cycles, the impulse intensity c need not be 
infinitcsimally weak provided that the inter impulse intervals are sufficiently long, in contrast to the conventional 
phase- reduction analysis of limit cycles driven by continuous signals [l(| HH, Q3 • It holds for largely fluctuating PRCs 
as well. By virtue of this fact, we can analyze desynchronization effect of random signals within the framework of a 
one-dimensional phase model, in contrast to Ref. [Tlj . as we discuss below. 



D. Linear stability of the synchronized state 

Whether synchronization occurs among the oscillators depends on the stability of the synchronized state. To 
investigate the linear stability of the synchronized state with the application of common impulses, we focus on the 
time evolution of a small phase difference tp ~ <j) — <f> between phase <j) an d </> of two nearby orbits. From Eq. ([8]) , the 
linearized evolution equation for tp is given by 

dtp(t) = [ tp-^-G((j},c)M(dt,dc). (12) 



Now we change variables to the logarithm of the absolute value of the phase difference, y = log| | . Using the stochastic 
chain rule for the jump process, we obtain 

dy(t) = ( [\og\^ + i)G\ct>,c)\-\ogmM{dt,dc) 
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= J \og\l + G'(<t>,c)\M{dt,dc), (13) 



where the prime(') denotes partial derivative by cp. The average growth rate of the small phase difference is charac- 
terized by the Lyapunov exponent A. In this case, it is defined as 



A = lim — los 



T->oc T 



= lim i f dy(t). (14) 



When A is negative, the initial phase difference decays exponentially so that the synchronized state is linearly stable. 
As usual, we postulate that the growing and shrinking of the phase difference is ergodic, namely, the long-time average 
slope of y(t) coincides with the ensemble average of its local slope, 

where (• • • ) denotes ensemble average over the marked Poisson process. Note that an individual increment dy(t) may 
exhibit discontinuous jumps of 0(1), but the ensemble average (dy(t)} is always of 0(dt). The expectation of the 
right-hand-side is calculated by replacing the dynamics of cf> with the single-oscillator stationary PDF p(<f) of <fi, as 



(dy(t)) = (J log\l + G'(<t>(t),c)\M(dt,dc) 

= [ dM4>) I log\l + G'(<j>,c)\(M(dt,dc)) 



= Xdt / d4p{<t>) / dcp(c) log \1 + G'(4>,c)\, (16) 



where in the second line, the ensemble average is separated into a conditional expectation with fixed <f> and average over 
p(cf>) because the integration variable and stochastic driving term are statistically independent. Thus, the Lyapunov 
exponent A is obtained as 

A = A / dc/> P (<t>) I dcp(c) log \1 + G'(ct>,c)\, (17) 



which generalizes the result obtained by random phase maps in Refs. [§] and [l3[ for general multiplicative coupling. 
The sign of the Lyapunov exponent depends on the shape of the PRC, G(<fi, c). When G'(<p, c) < —2 or G'((f>, c) > 0, 
the integrand, which gives the instantaneous growth rate of ip(t) a t is positive. Such regions tend to expand the 
phase difference between two orbits. When —2 < G'(<p,c) < 0, the integrand is negative, and the phase difference 
between two orbits tends to shrink. Linear stability of the synchronized state is determined by the overall balance 
between these two effects. 

For weak impulses, we can further simplify Eq. (|17[) by assuming that, on average, an oscillator is equally distributed 
on the limit-cycle, so p{4>) = 1 . This is a reasonable assumption in most cases where the effect of impulses are small [lU , 
as we discuss later in Appendix B. Under this approximation, the Lyapunov exponent A can be simplified as 

A = \[ d<t> [ dcp(c)\og\l + G'((j),c)\. (18) 

JO Jc 

Now, if the impulses are weak and the variation of the PRC G(<f>, c) is sufficiently small in such a way that — 1 < G'(4>, c) 
is always satisfied for all </> and c, i.e. when </> + G{<p, c) is a monotonically increasing function, A can be bounded 
from above as [HI, [l3| 

A < X f 1 dcf) f dcG'(4>, c)p(c) = A / p(c)dc [G{<j>, c)]g = 0, (19) 

Jo Jc Jc 

where we utilized the periodicity of G(cj>, c) in cf> and the inequality log(l + x) < x. Thus, for weak impulses, small 
perturbations are always statistically stabilized when averaged over the limit cycle, so that common impulses shrink 
the phase difference, irrespective of the details of the oscillator. A geometrical interpretation of this stabilization 
mechanism is given in Appendix C. 
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By a Taylor expansion of Eq. (|18|) , the Lyapunov exponent can be approximated for weak impulse as 
A = A J dcj) J dcp(c) 

where the first term drops out due to the periodicity of G(<j),c). At the lowest order approximation, G(<p,c) ~ 
Z(4>) ■ g(<f>, c) ~ J2k=i Zk(4 > )&k{4>> c )i the approximate Lyapunov exponent A is obtained as 

> At M r i 

A = -2 EE / wzttflzmMctt) 
k=i 1=1 ^° 

+2Z' k {4>)Zi{4>)M)M + z k {<\>)zM)K°'i)M)\ (21) 

for both Ito and Stratonovich pictures of Eq. |T]), where we introduced correlation functions (ofcO7) c (0) = 
J dcp(c)ak(4>, c)ai(4>, c), etc. As we derive in Appendix D, we obtain the same Lyapunov exponent in the diffu- 
sion limit of the jump process. 



G'(0,c)- 



G'(</>,c) 



J d<p Jjcp{c)G'(<t>,c) 2 < 0, 



(20) 



E. Clustered states 



If the PRC possesses a symmetry 

G(0,c)=G(<f> + ± 7 c) (22) 
m 

where m £ N and 1/m < 1 (i.e. m = 2, 3, 4, • • • ), the stability of the synchronized state (zero phase difference) also 
implies the existence of stable states separated in phase by 1/m. Defining tp as a small deviation from such a state, 
we set 

i/> = <f>' -</>-— . (23) 
m 

The linear stability analysis of a state separated in phase by 1/m yields the same Eq. (fT2|) when utilizing PRC 
symmetry, so that the resulting Lyapunov exponent A, Eq. (|18p . also takes the same value as that of the synchronized 
state. 

If A is negative, the phase difference between two orbits can stably take m different values. When many identical 
oscillators are driven by common impulses and also by small independent disturbances in such a situation, the 
oscillators will eventually split into m clusters. Any two oscillators inside the same group are synchronized, and 
those belonging to different clusters will take one of m — 1 phase differences, n/m where n = 1,- • • , m— 1. We call 
this an m-cluster state. As we demonstrate later, symmetry in the original limit cycle results in symmetry of the 
corresponding PRC, leading to cluster states. 



III. NUMERICAL SIMULATION 



In this section, we demonstrate synchronization, desynchronization, and clustering induced by common impulses by 
numerical simulations, and test our theoretical predictions quantitatively, using the FitzHugh-Nagumo (FHN) neural 
oscillator as an example. The effect of common random signals for uncoupled FHN oscillators have been discussed 
for Gaussian signals in Ref. [ll| and for a random telegraphic signal in Ref. [12j . However, quantitative comparison 
of the Lyapunov exponent predicted from the PRC with that directly measured from numerical simulations have not 
been fully done, especially in the desynchronization regime. In this section, we measure the PRC and the separation 
rate of nearby trajectories directly by numerical simulations, and quantitatively confirm the theoretical prediction 
based on the one-dimensional phase model. Results for the Stuart-Landau oscillator, which is qualitatively different 
from the FHN oscillator, is also given in Appendix E for comparison. 



A. Model 

For the simulation, we employ the FitzHugh-Nagumo (FHN) neural oscillator [28| as an example, described by the 
following set of equations: 



ii(t) = s(v + a ~ bu), 
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v(t) 



= v 



V 



3 



,3 



N(t) 

u + I + <j(v, c) 



(24) 



n=l 



Here, parameters e,a,b are fixed at e — 0.08, a = 0.7, 6 = 0.8, and we use the parameter Iq to control the oscillator 
characteristics. The last two terms of the equation for v describe impulses and noises. The h(t) is a common impulse 
with unit intensity, generated by a Poisson point process of constant-rate A. The function cr(v, c) describes w-dependent 
effect of the impulse to the oscillator (for simplicity, we do not consider the case where c takes multiple values in the 
following). The £(t) is a Gaussian white noise with intensity D describing independent disturbances to the oscillators. 

When both external disturbances are zero, a limit cycle exists for Iq 6 [0.331, 1.419], which is created by a subcritical 
Hopf bifurcation at either limits of Iq. Figure [T] displays portraits of asymptotic phase for two values of Iq in the 
absence of impulses and noises. At Iq = 0.8, the period of the limit cycle is T f=a 36.52, and the oscillator has a smooth 
phase portrait as shown in Fig. [IJa). Near the bifurcation point, Iq = 0.34, the period of the limit cycle isTw 46.79. 
The remnants of the destabilized fixed point exist at this parameter, as seen in Fig.QJb). We define the origin of phase 
cj) as the point where the variable v passes through the v = 0.9 line from below, where the FHN oscillator appears to 
emit a neuronal action potential. 



In simulation, we may use two algorithms corresponding to Ito and Stratonovich pictures of Eq. ([T]). If we treat 
the impulses as point events, allowing discontinuous jumps of the orbit of the oscillator, the results correspond to the 
Ito picture, Eq. (j4|). If we directly integrate short but non-zero- width continuous impulses, the results correspond to 
the Stratonovich picture, Eq. ©. When the effect of impulses arc multiplicative, the PRC is different between the 
two interpretations. Figure [3] displays the PRCs obtained using Ito and Stratonovich pictures for additive impulses 
(cr(w,c) = c) and for linear multiplicative impulses (ct(w,c) = cv). For additive impulses, both algorithms give the 
same PRCs. In contrast, for multiplicative impulses, it is seen that a difference in the pictures has a small but 
visible effect on the PRC. Of course, the PRC obtained using Wong-Zakai-Marcus approximation coincides with that 
obtained in the Stratonovich picture. All the following simulations are done using the Stratonovich picture of the 
original Eq. |T]), namely Eq. ((6]), using the Wong-Zakai-Marcus approximation of a continuous physical jump. 

At Jo = 0.8, the PRC for additive impulses is a smooth periodic function as shown in Fig. [Ha) for all values of the 
impulse intensity c, corresponding to the smooth phase portrait as shown in Fig.^Ja). In contrast, at Iq = 0.34, the 
PRC can become jagged as shown in Fig. [4Kb) for the impulse intensity c in a certain range. This reflects the existence 
of the unstable focus, which looks like a spiral in Fig. [ljb). If the impulse intensity is in such a range that the orbit 
on the limit cycle is kicked into this region, an initial phase difference can grow quickly because the asymptotic phase 
in that region varies so rapidly. 



Raster plots of N = 20 uncoupled FHN oscillators, Fig. [2 which indicate times that an oscillator passed through 
= 0, offer a qualitative picture of the phenomenon. Whether the impulse is introduced additively or multiplicatively, 
we find system and impulse parameters where phase synchronization occurs and where it does not. FHN oscillator 
has a large parameter range in Iq and c where common impulses cause the oscillators to synchronize, as shown 
in Fig. HJa). However, near the bifurcation, common impulses sometimes accelerate the desynchronization of the 
oscillators, as shown in Fig. [2Kb) . These distinct behavior can be understood by examining the PRC of the FHN 
oscillator at each parameter value. 

Figure [5] summarizes the relation between the shapes of the PRC and the dynamics of the oscillator ensemble in the 
phase space. For a smooth PRC with relatively small amplitudes obtained at Iq = 0.8 for additive impulse (tr(w, c) = c) 
as shown in Fig. EJa), the oscillators groups together on the limit cycle, namely, synchronize with each other, by the 
application of common impulses, Fig. [5Kb). For jagged PRC with strong amplitude fluctuations obtained by applying 
additive impulses near the bifurcation point Iq = 0.34 as shown in Fig. [5jc), the oscillators undergo desynchronization 
by the application of common impulses and scatter along the limit cycle, as shown in Fig.EJd). At the midway point 
between the subcritical bifurcation near Iq = 0.875, the limit cycle becomes symmetric about v = 0. If the impulse 
is applied in a linear multiplicative way {<j{v, c) = cv) in this Iq region, the PRC becomes doubly periodic, as shown 
in Fig. EKe), so that 2-cluster state appears as shown in Fig. [5]Jf). Even when the parameters of the oscillators are 
slightly inhomogeneous, these dynamical behavior remain qualitatively unchanged [44|. 

In the present case for FHN oscillators, synchronization gradually occurs whereas desynchronization occurs suddenly 
due to the narrow jagged part of the PRC, typically obtained near the bifurcation point of the dynamics. However, 



B. Phase response curves 



C. Synchronization, desynchronization, and clustering 
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we emphasize that desynchronization is not limited to such pathological situations. The PRC need not be rapidly 
fluctuating as long as intervals where the sufficiently steep slopes outweigh shallower slopes in Eq. (jTHJ) • As we see in 
the next section, our electrical oscillator is of this type, where desynchronization occurs even though the oscillator is 
far from a bifurcation and the PRC is smoothly. As an example of this from a well-known system, we include results 
from the Stuart-Landau oscillator in Appendix E. 

In Ref . , the mechanism for desynchronization is analyzed for FHN oscillators near the bifurcation point driven 
by finite strength white Gaussian forcing utilizing a two- variable phase-amplitude model to take into account the non- 
trivial transverse deviation from the limit cycle. In contrast, for our Poisson forcing case (and also for our previous 
case with random telegraphic forcing [l2| ) , we can eliminate the relaxation dynamics of the amplitude and isolate the 
effects of the phase perturbation, because the time-scale of oscillator relaxation back to the limit-cycle is assumed to 
be much shorter than the average impulse interval. We can thus stay within the framework of a single-variable phase 
model to describe both the synchronization and desynchronization, which elucidates the relation between the phase 
space structure (isochrons) and the desynchronization mechanism. 

This desynchronization mechanism is qualitatively similar to the mechanism of the singular behavior in circadian 
clocks proposed by Winfree [l8[ , who argued that the attenuation of the circadian rhythm by an external stimulus is 
due to the desynchronization of multiple independent circadian oscillators being kicked into the unstable focus of the 
limit cycle. In Ref. (l6j . application of such a desynchronization mechanism to neural populations are discussed. Based 
on the same idea, the sudden desynchronization seen in coupled chemical oscillators when a common perturbation is 
administered at the correct timing is also reported (l7j . Recently, Ukai et al. performed a very clear experiment 
of this phenomenon using genetically engineered photosensitive cells exhibiting circadian oscillations. 

D. Lyapunov exponents 

To quantitatively test our theory, we measured the Lyapunov exponent A in two ways, from numerically obtained 
PRC G(4>, c) using Eq. (fT8"|) and by observing the growth rate of the phase difference of two oscillators with period 
T using raster data, Fig. [2J for the case of common additive impulses. When the time difference At between the 
respective = events of two oscillators is below a threshold value (At < 0.02T), we converted the time difference 
into A</>(0) = At/T, and followed the evolution by measuring subsequent A(j>(t). This was done for many oscillator 
pairs, and for a given t, we calculated the average of log \Acf>(t)/ Acf>(0)\, as shown in Fig. [6] for Iq = 0.34. The slope 
of this line gives A. Figures [7] and [S] show that Lyapunov exponents measured from the simulation data show good 
agreement with those measured from PRCs as shown in Figs. HJa) and (b), respectively. 

E. On-off intermittency and switching 

While the instantaneous growth rate of the difference between two orbits, log 1 1 + G'(</>, c)|, may be a smoothly 
varying function with respect to phase </>, since the common impulses are received at random phases, the separation 
between two oscillators can be considered a random multiplicative process driven by the fluctuations in the growth 
rate around the average Lyapunov exponent A. In the absence of any disturbances, complete synchronization results 
if the average effect of a common impulse causes small phase perturbations to shrink as indicated by the average 
Lyapunov exponent over the limit cycle, Eq. (|18p . However, if small disturbances exist in the system, fluctuations in 
the instantaneous growth rate can occasionally amplify a small deviation in the system, causing intermittent transient 
desynchronization. 

Figures [9] (a), (c) show an example of the large excursions away from the synchronized and clustered state 
respectively, while Fig. EIb), (d) show the power-law PDF of laminar duration with the well-known exponent, 
— 1.5 0, Ho, HH, H3|. When the system exhibits clustering, the same mechanism leads to switching for sufficiently 
strong independent noises, with a power-law PDF of the life-times of the states as shown in Fig. [TO] 

This fluctuation is known as modulational or on-off intermittency, which was first discovered in a pioneering work by 
Fujisaka and Yamada [30l ] on a system of coupled chaotic oscillators, and later clarified to be an ubiquitous feature of 
many nonlinear dynamical systems with a symmetry. The intermittency typically arises in synchronization problems 
of dynamical elements, where a dynamical variable acts as a time-dependent fluctuating driving parameter for a 
second variable 0, [U, Ijffil . for example in the synchronization of chaotic lasers [33[ , in spin- wave instabilities [34| , and 
in nematic convection [35j . The same mechanism also applies to the present synchronization phenomenon induced by 
common random signals, where the phase difference is multiplicativcly modulated by rapid random forcing due to the 
common random signals. The statistical properties of the separation of trajectories in such a situation was previously 
investigated in detail in Ref. [36[ for two uncoupled maps receiving a common noisy drive, a system very similar to 
the one currently under consideration. 
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IV. EXPERIMENT 

In this section, we present the results of our experiments on an electrical limit-cycle oscillator. S ync hronization of 
electrical limit-cycle oscillators induced by a continuous random signal has been realized e.g. in 13811 . but without 
quantitative comparison with the theory. The deduction of the PRC of noisy neural oscillators [39[ and the use 
of the PRC in predicting oscillation stability, or firing reliability, for cells receiving complex stochastic input (ioj 
have also been discussed in the neuroscience literature recently, which mainly focus on the synchronizing effects of 
common random signals. In this section, we experimentally measure the PRC in our electric circuit experiments, and 
quantitatively compare the Lyapunov exponent theoretically predicted from the PRC with that measured directly 
from experimental time sequences, in both synchronization and desynchronization regimes. 

A. Setup 

To observe common impulse-induced synchronization and desynchronization, we experimented on an electrical limit- 
cycle oscillator. Figure [TT] shows the circuit diagram and limit cycle for our experimental system, a battery-powered 
LED-flashing oscillator, where the natural period of the oscillator is T w 0.79s. Voltages were measured at two 
locations in the circuit, Chi and C/i2- These voltages were taken to be phase space variables with which we measured 
the limit cycle. The limit cycle displayed slight wandering in the phase space of V c hi x Vchi as the natural frequency 
of the oscillator drifted on the order of 1% over the course of the experiment, so a new = line in phase space was 
uniquely chosen every 5 experimental trials after which the limit cycle was re-calibrated. This drift in frequency is 
equivalent to a small non-identicality of oscillators in an ensemble. The = line was drawn perpendicular to the 
region of the limit cycle where the oscillator displayed the fastest dynamics, as such a choice ensured that the = 
crossing event could be measured with the least uncertainty. 

The impulses were created by computer, which delivered a voltage signal V g via the output of a data acquisition card 
to the gate of the metal?oxide?semiconductor field-effect transistor (MOSFET) Mi acting as a constant current source, 
i.e. a state- independent, additive impulse. In order to simulate an ensemble of identical oscillators receiving common 
impulses, the experiment was repeated many times employing an identical train of impulses throughout the trials 
with either random or identical initial conditions to investigate synchronization or desynchronization, respectively. 

B. Phase response curves 

As in the simulation, we obtained the PRC experimentally by applying impulsive stimulus to the circuit. We varied 
the PRC characteristics by varying the location of the circuit to which impulses were applied. Figure [T2] shows the 
experimental PRCs Gi(0, V g ) and G2(0, V g ) obtained by stimulating Chi and Ch2 for several impulse intensities. 
Qualitatively different PRCs were obtained when we stimulated different locations. In contrast, varying the stimulus 
intensity to the same location resulted in similar but systematically different PRCs. 

An experimental measurement has many sources of noise, which seriously degrades the first derivative calculated 
from a discretely sampled PRC, which is necessary in calculating the Lyapunov exponent. Assuming that the fluctua- 
tions seen in the derivative of the PRC arc due to experimental noise and that the underlying response of the system 
is smoothly changing with respect to phase, we must infer the underlying smooth response from our noisy data. 

To this end, we smoothed the PRC using a non-Gaussian Kalman filter [37| based on a Markov state space model 
with a transition probability distribution given by the Cauchy distribution. The Kalman filter iteratively finds the 
"true" values, {0i, 02, ■ • • , <Pn}, at each i of the data given the observed data {0i, 02, • • • , 0jv} such that conditional 
probability 

p{{(t>i}\{4>i} , {filter parameters}), (25) 

is maximized. The method may be combined with Bayesian methods such as the Expectation-Maximization algorithm 
to choose the best filter parameters. We did not perform this step, but chose parameters that preserved the global 
shape of the PRC while yielding a smooth first derivative. This smoothed PRC was then used in Eq. (fT5|) to calculate 
the Lyapunov exponent. 

C. Synchronization and desynchronization 

Figure 02] shows the voltage traces at Chi for several different trials showing electrical oscillators synchronizing 
and desynchronizing due to common impulses. The PRC obtained when impulses were applied to Chi, Fig. II2f a). 
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had a large jump, and the effect of the common impulse was predominantly a dcsynchronizing one. We were able 
to detect a synchronization-desynchronization transition in Ch\ by varying the impulse intensity (data not shown), 
but because the synchronizing effects were relatively weak, we chose instead to investigate the synchronization due to 
impulses applied to C/12, where the PRC is smooth as shown in Fig. [T2Tb) so that we were able to observe common- 
impulse induced synchronization clearly, as shown in Fig. I14f a). We verified the theoretical predictions by measuring 
the Lyapunov exponent using the two methods outlined in the previous section. The results of the two methods of 
measurement of the Lyapunov exponent arc summarized in the caption of Fig. 1151 Considering the frequency drift 
among trials, the agreement between the results of the two methods seems reasonable. (45j . 

Note that our electrical limit-cycle oscillator is not near the bifurcation point and the PRC is not jagged in contrast 
to the FHN oscillator, even though the oscillators are desynchronized when impulses are applied. Therefore, the 
desynchronization effect is rather gradual, as shown in the raster plot, Fig. I14f b). which is qualitatively similar to the 
situation with strong impulses on the Stuart-Landau oscillators that we discuss in Appendix E. 



Through theoretical analysis, numerical simulations, and circuit experiments, we have demonstrated that phase 
synchronization, desynchronization, and clustering can be realized in a system of identical, uncoupled oscillators 
acted on by common Poisson impulses by changing the parameters of each individual oscillator, and also by varying 
the strength of the common impulse. We have clarified that once the phase response curves of the oscillator is known, 
such phase coherence phenomena can be quantitatively understood in a unified way. The desynchronization is an effect 
of the finite size of impulses that we use, in sharp contrast to the exclusively synchronizing effects that infinitesimal 
perturbations have. 

It is quite remarkable that the synchronization and desynchronization of an ensemble of oscillators may be controlled 
in this way through a random signal. Depending on the oscillator characteristics, we could alter the coherence 
properties of the ensemble simply by changing the intensity of external random impulses. This approach would have 
potential applications in which it is desirable to change the global coherence properties of a network of oscillators. 
Applying common impulse of suitable intensity, it would be possible to effectively "switch off' or "switch on" the 
synchronization without changing the oscillator characteristics or modifying the coupling strength. 

In this paper, we characterized the local stability of the synchronized states by a Lyapunov analysis, but not the 
global stability. In Ref. [l4| , we have presented a global stability analysis of the system for common Gaussian- white 
driving using an averaging technique for nonlinear oscillators. Similar formulation is also possible for the present 
common impulsive driving. The clustering phenomenon was realized only numerically for the FHN system, because 
our present experimental circuit does not appear to have a symmetric limit cycle or PRC compatible with clustering. 
Other limit-cycle electric oscillators with nearly symmetric limit cycles do exist, and investigation of clustering with 
such systems is now under progress. We expect that on-off intcrmittency and switching can also be observed in such 
electric oscillators. Detailed analysis on these topics will be reported in the future. 



The authors wish to thank Y. Kuramoto for helpful discussions, especially suggesting us to analyze the impulse- 
driven case, and H. Fujisaka for all his insightful comments and advice. We also thank K. Aihara, A. Uchida, K. 
Yoshimura, H. Suetani, T. Shimokawa, J. Teramae, D. Tanaka, T. Kobayashi, Y. Kawamura and Y. Tsubo for 
various useful comments and information. This work is supported partially by the Grant-in- Aid for the 21st Century 
COE "Center for Diversity and Universality in Physics", and partially by the Grant-in- Aid for Young Scientists (B), 
19760253, 2007, from the Ministry of Education, Culture, Sports, Science, and Technology of Japan. 



In this paper, we model the process of an oscillator receiving impulses as a Poisson-driven Markov process, or jump 
process [22|, |23j. A general stochastic process X(t) driven by Poisson random impulses, 
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APPENDIX A: POISSON-DRIVEN MARKOV PROCESS (JUMP PROCESS) 



N(t) 




(Al) 



71=1 
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with point-like impulses h(t) should be interpreted in the Ito picture as 

N(t) 



X(t) = X(0)+ f F(X(s))ds + VG(X(t„-0),c n ), (A2) 
J n=l 

which is described as an integral equation of the form 

X(i) = X(0)+ [ F(X(s))ds+ [ [ G(X(s),c)M(dt,dc), (A3) 

JO JO Jc 

or as a jump stochastic differential equation (SDE), 

dX(t) = F(X)dt + J G(X, c)M(dt, dc). (A4) 

Here, M(dt,dc) represents a Poisson random measure [22I , [23l |, which gives the number of incident points during 
[t, t + dt] having mark (jump magnitude) in [c, c + dc], and N(t) = J Q J c M(dt, dc) is the number of incident points 
during [0, t]. As usual in Ito-type stochastic differential equations, functions of X(t) and the Poisson random measure 
M(dt,dc) at the same instant of time are independent. The expectation of M(dt,dc) is 

(M(dt, dc)) = \dtp(c)dc, (A5) 

where A is the rate of the Poisson process and the marks are distributed as p(c). Similarly, the covariance of M(dt\, dc\) 
and M(dt 2 ,dc 2 ) is 

Cov[M(dti,dci),M(dt 2 ,dc 2 )} = \5(t 2 - t 1 )dt 1 dt 2 p(ci)5(ci - c 2 )dc 1 dc 2 . (A6) 

As in the case of Wiener-driven Ito stochastic differential equations, we need to use a special rule in calculating the 
differential of a transformed process. A transformed process V(t) = V(X(t)) obeys a jump SDE of the form 

dV{t) = [grad x y(X) ■ F(X)} dt + J [V(X + G(X, cj) - V(X)]M(dt, dc), (A7) 

which is a stochastic chain rule for the jump process. The chain rule changes an additive process into a multiplicative 
process, just as in the case of white- Gaussian driven Markov processes. 

APPENDIX B: SINGLE OSCILLATOR PHASE DISTRIBUTION 

In the calculation of the Lyapunov exponent, we simplified the calculation by assuming that a single oscillator 
receiving random Poisson impulses is evenly distributed in phase when averaged over a long period of time. This is 
strictly not the case, as can be seen from the map (f)(6) = 9 + G(6,c). It is obvious that certain values of (f)(9) are 
arrived at more often than other values, because G(9, c) is generally a nonlinear function of 9. For completeness, 
we calculate the deviation from a flat phase PDF by the forward Kolmogorov equation [23[. To this end, we define 
r)((f>, c) = G(9, c), which is the jump as a function of the destination coordinate <j>. Then the probability density p((f), t) 
of </> obeys 



d d 

—p(4>, t) = -Qi [up(<t>, t)\ - Ap(& *) + A 



dc. (Bl) 



We find the stationary PDF by setting dp((f>, t)/dt = 0. It is clear that the combination e = X/ui determines the effect 
of the impulses, which is a small parameter from our initial assumption (the rate of the Poisson process A is small). 
We thus expand p(4>) in powers of e away from the stationary PDF as p((f>) = 1 + epi((f>) + • ■ • • Substituting this into 
the forward Kolmogorov equation, the first non- vanishing terms are of order 0(e 1 ): 



= ^ 



Pi((f>) + 1- / p( 



dc. (B2) 



If the PRC does not change too rapidly, —1 < dri((f>,c)/d(f) < 1, and the absolute value of the integrand may be 
removed, yielding 

§jPM = - jp(c)^n(^ c)dc. (B3) 
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Then 

pi(0)=C- [ p(c)ri{<f>,c)dc. (B4) 

J c 

Since the integral of p(<fi) over one period of the 0(e ) term is 1 due to normalization, higher order terms 
Pi (0)jP2(0)i ' ' ' must vanish upon integration over a full period. Specifically, the J pi(<f>)d<f> — condition yields 



where 



C = j p{c)rj(c)dc (B5) 

J c 



77(c) = / v (<j>, c)d<f>. (B6) 







Therefore, we have for the first order approximation to the stationary phase PDF, 

p{4>) = 1 + e J [77(c) - n{<f>, c)] p(c)dc + 0(e 2 ). (B7) 

As we see from Fig. [TT1 p(</>) is close to constant as long as the parameter e is small. In Ref. [l3[, we argued that 
for Poisson impulses, the lowest order correction to the uniform phase density does not contribute to the Lyapunov 
exponent based on a perturbation expansion of a Frobcnius-Perron-type equation. Here we only point out that the 
correction to the uniform density is of 0(A/u>), which is small if the Poisson rate A is sufficiently smaller than the 
oscillator natural frequency uj. 



APPENDIX C: GEOMETRIC INTERPRETATION OF THE SYNCHRONIZING MECHANISM 

Weak impulses always synchronize uncoupled oscillators. Here we show through a simple Floquet analysis that this 
synchronization is a consequence of the stability of the limit cycle against weak perturbations. The generic oscillator 
under consideration in this article is described by an ordinary differential equation of the form 

X(t) = F(X), (CI) 

with a stable limit-cycle solution Xa(t) with period T. Linearizing Eq. (|C1[) with respect to a small perturbation u(t) 
from the limit cycle, we obtain 

u{t) = DF(X)\ x=Xa{t) u (C2) 

where DF(X)\ x =x (t) is a periodic MxM Jacobian matrix. Ordinary differential equations of this form with periodic 
coefficients have solutions of the form u(t) = Q(t)e Rt u(0) where R is a constant MxM matrix, and Q(t) = Q(t + T) 
is a T-pcriodic MxM matrix. Since Q is periodic, we have u(T) = e RT u(0). There is a corresponding value of 
the constant matrix R for each point on the limit cycle. The eigenvalues and eigenvectors of e RT , {Xi} and {e^} 
respectively, have the property that e\ is in the direction along the limit cycle, Ai = 1 and |A,| < 1 for i £ {2, • • • , M}. 
The eigenvectors are known as the Floquet eigenvectors, and each point on the limit cycle has its own set of Floquet 
eigenvalues and eigenvectors. 

Figure [TBI shows two infinitesimally separated orbits, 1 and 2, on the limit cycle at a, Xo (</>), and b, Xo(<f>) + z(0), 
with a phase <p isochron passing through a. Now let the two oscillators receive a common additive impulse c at t = 0. 
The two oscillators jump discontinuously to a at Xq{4>) + c and b at Xo(<f>) + z(0) + c, with a phase </> isochron passes 
through a. The set of Floquet eigenvalues and eigenvectors at <j) an d <j> are ({A^}, {e^}) and ({A^}, {e.;}), respectively. 

Since the impulse is additive, ab=ab= z(0). Expanding the difference vector z(0) by the Floquet vectors, we obtain 
z(0) = aiei = J2i OLi&i- It is then obvious that |ai| > |6i|. 

If the oscillator continues unperturbed for one period, the oscillator that received the impulse at a at t = will now 

be at a*, and z(T) = e RT z(0) = J^i Ka^i- Since Ai < 1 for i e {2, • ■ ■ , M} from the stability of limit cycles, all 

components of z(T) with i > 2 will shrink, leaving only the first component along the limit cycle. We thus see that 
|z(T)| < |z(0)| by virtue of |ai| < |ai|, namely that the application of a common impulse always shrinks the small 
separation between two orbits. 
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APPENDIX D: DIFFUSION LIMIT 



We here derive the diffusion (Gaussian-white) limit of the impulse-driven oscillators for weak and frequent impulses. 
In taking the diffusion limit, we see that common impulse always results in the synchronization of oscillators since the 
Lyapunov exponent is bounded above by 0, and we see that desynchronization can only be understood by analyzing 
finite-magnitude perturbations to the limit cycle orbit. 



1. Diffusion limit 



We have until now considered finite Poisson impulses whose inter-impulse times are much longer than the natural 
period of the oscillators. The condition for the phase reduction is also satisfied by setting the combination of impulse 
intensity and the Poisson rate appropriately small, so that we may also consider a situation where the effect of Poisson 
impulses become infinitesimal but with inter-impulse times that are much faster than the natural time scales of the 
oscillators. Specifically, we consider the limit A — > oo and the effect of impulses <7fc — > such that \((Jk&i)c is kept 
constant and higher order terms like \(crkcricr m ) vanish (k,l,m = 1, ■ • ■ , M). To take this limit, it is necessary that 
the net effect of the impulses vanishes, namely, 

(o*(0,c))c = J p(c)a k (X (ct>),c)dc = 0, (fc = l,--- ,M), (Dl) 

where we introduced the notation (A(<f>, c)) c = J A(c/>, c)p(c)dc with fixed <fi. Under these conditions, we can consider 
the diffusion limit of a stochastic jump process. 

The Kramers-Moyal expansion of the Chapman-Kolmogorov equation is [4l| 



d P (<p, t ) _ ^ i f d 
dt 



^.{-96 K^M^j, (D2) 



ra=l 

where the Kramers-Moyal coefficient is given by 

y ' At->o At v ; 

Here, A(f> is the jump within duration At of the stochastic process starting from cp, and the conditional average is 
taken over possible realizations of the stochastic process starting from (j>. If the coefficients higher than the second 
order vanish, we obtain a Fokker-Planck equation 

= ~lk M^'*)] + \w W)p ^ t)] (D4) 

for the Wiener-driven Markov process, whose drift v{4>) an d diffusion coefficient D(<f>) are given by 

v(4>)=KW (</>), D(4)=KW^). (D5) 

We now find (A0) and (A</> 2 ) for a jump process, where the expectation is to be taken with fixed <f>. Starting with 
Eq. ®, 

At 

2^ 



Similarly, 



(A^) = uAt+ G(4>,c)(M(dt,dc)) +0(At 2 ) 

JO Jc 

= LuAt + \At(G((f>,c)) c + 0(Ai 2 ). (D6) 
(Acf) = ((A0-(A0)) 2 )+O(At 2 ) 

"At pAt r- r 

/ / / G(<t),c)G((t)' ,c') Cov[M(dt,c),M(dt',c')} + 0(At 2 ) 

JO Jc Jc' 

XAt(G(4>,c) 2 ) c + 0(At 2 ). (D7) 
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It can be checked that (A<fi 3 ) and higher-order moments vanish by taking the diffusion limit, so that they can be 
dropped. The Kramers-Moyal coefficients are given by 



tf«(0)=u, + <G(^c)> Cj A^>(0) = (GW>,c) 2 ) c , *f<»>(0)=O (n>3). 



(D8) 



Hereafter, for simplicity, we may not explicitly indicate the dependence of the function g(Xo(<p),c), <r(Xo((p)), or 
Z((j>) on </> and c. To evaluate the v(</>) and D((f>), we must evaluate (G(</>, c)) c and (G((j),c) 2 ) c . Since we assume the 
effect of impulses g to be small, we first rewrite G(</>, c) by Taylor expanding, 



G(cf>,c) = 0(X o (0+ff(Xo(0,c))-0 

M M „ 2 

5EE 



M 

E 

fe=i 



2 ^ — ' ^ — ' dXhdXi 

X=X {<j>) fe=1 ;_i K 1 



3fe3/ 



X=X o (0) 



^From the definition of the phase sensitivity function, Eq. (jlip . we obtain 

d(j> 



dX k 



Zk(<f>), 



x=x {<f>) 



and 



d 2 



dX k dXi 



dZ l 



dX k 



dZi dej) 



X=X o (0) 



(D9) 



(D10) 



(Dll) 



X=X o (0) 



X=X (» 

Keeping only terms up to second order in g, we obtain 

(G(0,c)) c = Y. z ^9k) c + \Y. z ^km) c , 



(G(0,c) 2 ) c = Y, Z k z > 



l\9k9l) 



(D12) 



A./ 



Depending on the picture of the original Eq. ([I]), the approximate jump magnitude g(Xo(0), c) for a given mark c is 

<t(X o (0),c) (Ito), 
<7(X o (0),c) = { . A . (D13) 



IX 



x=x 



(Stratonovich). 



For v(4>) and -D(</>) in the Ito picture of Eq. (p}, we obtain 

v{<j>) = u J + ^Y.ZkZ' l {a k 



k.l 



D{4>) = Y,ZkZi{a k cji) 



(D14) 



k.l 



where we utilized the assumption {a k {4>,c)) c = J c p(c)a k (Xo((j)),c)dc = 0. For v{4>) and D{4>) corresponding to the 
Stratonovich picture of Eq. ([1]), we obtain up to 0(a k ai), 



1 



1 d(j 

g k {Xo((/)),c) DX k \ x =x (<f>) + ^D 1 X k \ x =x Q (<t>) = <?k + jE^Sf ' 



(D15) 



g k (X ((j)),c)gi(Xo(4>),c) = a k a u 
so in the Stratonovich picture of Eq. ([1]) , we obtain 

A 



(D16) 



v(<f>) 



^ {Z h Z[{<j k (ji) c + Z k Zi{a' k ai) c ) , 



A.; 
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D {4>) — ^ Z k Zi(a k (Ji) c , 



where we used 



k.l 



dak _ dak dej) 
ai dXi ~ ° l d(j> dXi 



dak 



(D17) 



(D18) 



Now that we have v(4>) and D{<j>) for both cases, we write an FPE, and find the corresponding Ito SDE. The Ito SDE 
corresponding to the Ito picture, Eq. (0J), reads 



<ty(t) = iw + - ]T Z k Z[{a k a x ) c dt + Va]T Z k G km dW m (t), 



(D19) 



kj 



k.m 



where we introduced M independent Wiener processes {dW m (t)} (m = 1, • • ■ , M) |4a | and an M x M coupling matrix 
Gk m {4>) that satisfies 



^ Gkm(<f>l)Gl m (<h) = (cfct^lW'^c- 

m 

Similarly, the Ito SDE corresponding to the Stratonovich picture, Eq.®, leads to 

A 



(D20) 



d<j)(t) = 



- {Z k Z[{(T k (Ti) c + Z k Zi{a' k ai) c ) 



k.i 



dt + ^\J2z k G km dW m (t). 



(D21) 



Using the transformation rule between Ito SDE and Stratonovich SDE [26|,|27|, this equation can concisely be expressed 
as a Stratonovich SDE 



(S) = udt + Vx^2 z k (4>)G km (4>)dW m (t), 



(D22) 



which was the starting point of the previous works [l(J Q3 ■ 



2. Linear stability of the synchronized state 
In both pictures of Eq. (fTJ), the diffusion-limit Ito SDE takes the form 

d(f>(t) = [u + a{4>)\ dt + J2 bm(<f>)dW m (t), 



(D23) 



where a(<f>) is periodic in <p, and b m (4>) = Zk(4>)G kr n{4>)- We are interested in the linearized dynamics of the 

small perturbation tp(t) to <f>(t), 



dip(t) 



(D24) 



Using the Ito formula [26|, (27| for changing variables to y = log \ tp\, 



dy(t) 



dt + J2Kni<P)dW m (t). 



The expectation is calculated by replacing the dynamics with the single-oscillator phase PDF, p(<f) 
Lyapunov exponent is given as 



(D25) 
1, so the 

(D26) 
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where the integral of a'(<j>) vanishes due to the periodicity of a (0), and the noise term vanishes because b' m ((j)) and 
dW m (t) are independent in the Ito SDE and the expectation of dW m (t) is 0. Therefore, the Lyapunov exponent is the 
same no matter the picture of the SDE Eq. ||7J). Inserting &,„(</>) = vA^fe Z k (4>)Gkm{.<l>), the summation in Eq. (|D26|) 
can be calculated as 



E & " = x J2 z 'k 



k.i 



GkmGlr, 



k,l 



Zi 



k.l 



k.l 



where we used 



etc., so we finally obtain 



-^E (Z' k (p k oi) c Z[ + Z' k {a k< rl) e Zi + Z k {a' k a t ) C Z[ + Z k {p' k a\) B Z l ) , 



d0i 



A = ~ E / \ Z 'M°W)* + M' k ZlM)e + Z k Z X (p' k o\) c \ 
k.l J ° 



(D27) 



(D28) 



(D29) 



This expression coincides with the approximate Lyapunov exponent that we obtained by a Taylor expansion in Eg . (I20[) , 
and gives a multiplicative generalization to the previous results obtained by Teramae and Tanaka in Ref. [l0( (our 
result in Ref. [l4| includes this result). 



APPENDIX E: STUART-LANDAU OSCILLATOR 



In Sec. Ill, we discussed the case in which the response of the oscillator to sufficiently strong perturbations result 
in PRCs that appear jagged, Fig. [4Jd). For such PRCs, the desynchronization is intuitive: if a nearly-synchronized 
group of oscillators near such a jagged response receive an common impulse, they end up with widely distributed 
phases [lU, [l3| . In Ref. [ll[ , the same situation is described differently, where the importance of the " heavy-tails" of 
the distribution of relaxation rates of transverse perturbations for oscillators near the bifurcation point is emphasized. 

However, the PRC need not have such a pathologic shape for desynchronization. It may even be sinusoidal as shown 
in Ref. Sec. 15. Such a case occurs with the Stuart-Landau (SL) oscillator, which describes the small-amplitude 
oscillations near the supercritical Hopf bifurcation point of a general system of ODEs [2lJ ■ 

Consider the following SL oscillator driven by random Poisson impulses: 

N(t) 

u = (u — cqv) — (u — C2v)(u 2 + v 2 ) + a(v, c) h(t — t n ) 

n=l 

v = (v + cqu) - (v + c 2 u){u 2 + v 2 ), (El) 

where h(t) and a(v, c) as described above for the FHN oscillator. For comparison with the FHN oscillator, we 
followed the same procedure for the SL oscillators and found the existence of synchronizing and desynchronizing 
impulse strengths (raster data not shown but are qualitatively similar to Fig. [2]). PRCs are shown in Fig. [TWa) 
and Lyapunov exponents obtained using Eq. (|18p and by direct measurement are shown in Fig. I18f b). The PRCs 
are almost sinusoidal but slightly deformed because the impulse intensity c is finite. As expected, the Lyapunov 
exponent A(c) shown in Fig. fT8T b) is qualitatively very similar to that obtained in Ref. Q calculated for a circle map 
with a sinusoidal PRC receiving common Poisson impulses, which predicts synchronization for weak impulses and 
desynchronization for stronger impulses. 

Note that in our present treatment of Poisson-driven limit cycles, we do not need to discuss the FHN-type oscil- 
lator and the SL-type oscillator separately. We can simply adopt the same one-dimensional phase model with the 
standard definition of the PRC, which quantitatively predicts the Lyapunov exponent in both synchronization and 
desynchronization regimes. 
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FIG. 1: (Color online) Asymptotic phase for FHN oscillator with limit cycle in black. To = 0.8 and 0.34 in (a), (b), respectively. 
The center of the spiral in (b) occurs at the intersection of the nullclines, and is the remnant of a destabilized fixed point as 
the oscillator passes through a subcritical Hopf bifurcation where Iq is the bifurcation parameter. 
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FIG. 2: Raster plots of FHN oscillators showing synchronization and desynchronization due to common impulses with rate 
A = 1/4T for both. Time axis normalized by natural frequency of oscillators. Each tick indicates the time oscillator passed 
through <p — near the limit cycle, (a) Synchronization of FHN (Jo = 0.8, c = 0.3, D — 5 x 10 -6 ) and (b) Desynchronization 
of FHN (Jo = 0.34, c = 0.2, D = 5 x 10 -8 ). The wide blank without ticks in (b) corresponds to the situation where the orbit 
is transiently trapped around the unstable focus. 
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FIG. 3: (Color online) Comparison of Ito vs. Stratonovich interpretations of Eq. JT} on the PRC G((j>, c) of FHN for an impulse 
whose jump size is c = 0.2. (a) Additive impulse (a(v,c) = c) and (b) Linear multiplicative impulse (a(v,c) = cv). The curve 
"Ito" is calculated by affecting a discontinuous jump, i.e. impulse duration is 0. The curve "Stratonovich" and "Strat." is 
calculated by continuously changing v using a narrow rectangular waveform of temporal width 0.0002. The curve "Exp." is 
calculated using the Wang-Zakai-Marcus approximation for the continuous narrow impulse, namely, discontinuously changing 
the orbit by an amount g(v,c) — (e c — l)v. 
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FIG. 4: (Color online) (a) PRCs for FHN with 7 = 0.8 and additive impulse intensities c G [-0.4, 0.4], (b)PRCs for FHN with 
Io = 0.34 and additive impulse intensities c = —0.20, 0.03, 0.20. 
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FIG. 5: (Color online) Various PRCs and coherent states of FHN oscillators, (a), (c), (e) show PRCs G(<j),c), and (b), (d), 
(f) show corresponding phase-space diagrams of 200 oscillators a sufficient time after the initial condition shown in the insets. 
Poisson rate is A = 1/4T for (a), (b), and (f), and independent noise is D = 5 x f0~ 6 for (b) and (f), and D = 9 x f0~ 9 for (d). 
The control parameter is To = 0.875 for (a), (b), (e), (f), and Jo = 0.34 for (c) and (d). For weak additive impulse, the PRC 
G{4>,c) is a periodic function, (a), and the 1-cluster (synchronized) state, (b), appears. At Jo = 0.34, the PRC G(<j>,c) becomes 
jagged when the additive impulse intensity is in a certain range as shown in (c), which often leads to common- impulse induced 
desynchronization, (d). For multiplicative impulse at Jo = 0.875, a doubly periodic PRC G(<f>,c), (e), leads to the 2-cluster 
state, (f). 
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FIG. 6: (Color online) Pairwise growth times of perturbations to FHN at Io = 0.34 for several values of the impulse intensity 
(A: c = 0.03, B: c = -0.6, C: c = 0.4, D: c = 0.2, E: c = -0.35, and F: c = -0.2). Data was taken for 100-200 trials, each 
with an ensemble of 20 oscillators, and with Poisson impulse rate A = 1/4T. Slope of linear least-square fit gives the Lyapunov 
exponent. 
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FIG. 7: (Color online) Comparison of the Lyapunov exponents A between the direct measurement from the raster plot and the 
theoretical prediction from the PRC for FHN with parameter Jo = 0.8 driven by additive impulses of intensity c, and Poisson 
impulse rate A = 1 /AT. 
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FIG. 8: (Color online) Comparison of the Lyapunov exponents A between the direct measurement from the raster plot and the 
theoretical prediction from the PRC for FHN with parameter 7o = 0.34 and additive impulses with rate A = 1/4T. Labels A, 
B, • • ■ , F correspond to those in Fig. [6] 



27 



a) 



0.10 



<f 0.05 
0.00 



¥ 








i 









2x10 




2x10 



c) 
P 



d) 
P 



10~ 

io 2 
io ] 

10 



o 



- 1 — L L I.I LIU 1 1 .1 J 1.1 IJ - 

' -NJ -1.5 


• 

1 — r r ri nr 


i — i »mwr 



io 1 io 2 io 3 

t/T 



10~ 

io 2 
io 1 
io c 



j 1 — L L I.I Lib 


i — i i j j i.i u . 

-1.5 : 


• 

- j 


!*\ ■■; 

• «w* mm - 



io 1 io 2 io 3 



t/T 



FIG. 9: (Color online) On-off intermittency exhibited by 2 oscillators in synchronized (Jo = 0.8, additive impulses with 
c = 0.1), and clustered (Io = 0.875, linear multiplicative impulses with c = 0.1) states. Independent noise with D — 9 x 10~ 9 , 
and impulses with A = 1/4T were used. Phase difference A(/> between the oscillators from stable configuration is small (laminar 
region) much of the time, but large occasional bursts occur, (a) long-time evolution of A<j>(t), which shows excursions away 
from the synchronized state, (b) distribution of laminar duration corresponding to (a) (arbitrary normalization), (c) long-time 
evolution of A<j>(t) from the 1/2-out-of-phase clustered state, and (d) distribution of laminar duration corresponding to (c) 
(arbitrary normalization). Oscillators are considered to be in the laminar state when A(f> < 0.0013 away from synchronized or 
clustered states. Laminar distributions exhibit power laws with exponent —1.5. At this weak independent noise intensity, the 
phase difference between the oscillators takes only either or 1/2 depending on the initial condition, and switching between 
the clustered state occurs is a very rare event. 
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FIG. 10: (Color online) Transition between clustered states for 2 oscillators. A(/> is the phase difference between the two 
oscillators. Poisson impulses with rate A = 1/4T and c = 0.1 was used. A larger independent noise with D — 3 x 10 -4 is added 
in order to facilitate the transitions between single and 2-cluster states. 




FIG. 11: (a) Diagram of electrical circuit with limit-cycle behavior. Computer-generated impulses control V g , which turns 
MOSFET Mi current source on/off. Switch Si allows us to send common impulse to either Chi or C/12. (b) Limit cycle of 
electrical circuit produce by measuring voltages at Chi or C/12 as given in (a). 
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FIG. 12: (Color online) PRCs Gi(<f>,V g ) and G2{4>, V g ) of electrical oscillator obtained by stimulating (a) Chi, and (b) Chi, 
which show responses of oscillators that desynchronize and synchronize, respectively, upon receiving common impulses. Each 
curve is labeled with a letter (A, B, C, ■ ■ ■ ) that corresponds to a location where impulse was applied (Chi or Chi), and the 
MOSFET gate voltage creating the impulse, which corresponds to the Poisson mark c. 
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FIG. 13: (Color online) Representative waveforms of electrical oscillators undergoing common-impulse induced synchronization 
(V g — —7.79V added to Ch-i) (a) and desynchronization (V g = —7.69V added to Chi) (b) measured at Chi. The Poisson 
impulse rate is A = 1/4T. Voltages traces in (b) actually extend below 0.2V, but have been clipped to show detail. 
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FIG. 14: Raster plots of electrical oscillators showing synchronization and desynchronization. The Poisson impulse rate is 
A = 1/4T. Each tick indicates the time oscillator passed through <f> = 0. (a) Synchronization of electrical oscillators (impulse 
added to C/12, V g = — 7.79V) and (b) Desynchronization of electrical oscillators (impulse added to Chi, V g = —7.64V). 
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FIG. 15: (Color online) Growth times of perturbations to electrical oscillator, where (A, B, C, • • • , F) correspond to that 
introduced in Fig. 1121 Data was taken for 80-150 trials, each with an ensemble of 20 oscillators, with Poisson impulse rate 
A = 1/4T. A comparison of the Lyapunov exponent with theory is shown in the table. 
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FIG. 16: (Color online) Schematic of evolution two nearby orbits, a = Xo(4>) and b = Xo(4>) + z(0) represent the spatial 
points at which 2 orbits receive a common additive impulse c on the limit cycle. The oscillators jump to a = Xo(<j)) + c 
and b = Xo(4>) + z(0) + c. After the oscillator completes one period of unperturbed motion, through analysis of the Floquet 
eigenvectors and eigenvalues, we see that the difference vector has shrunk, i.e. |«(r)| < |z(0)|. 
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FIG. 17: (Color online) The stationary phase distribution of a FitzHugh-Nagumo oscillator receiving random Poisson impulses, 
Jo = 0.8, c = 0.5 and A = 1/4T, calculated using a perturbation expansion of the forward Kolmogorov equation and by direct 
numerical simulation. 
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FIG. 18: (Color online) a) PRCs for SL model (co = — C2 = 12) for various additive impulse intensities, b) Comparison of the 
Lyapunov exponents A between the direct measurement from the raster plot and the theoretical prediction from the PRC for 
SL oscillators driven by additive impulses of intensity c, and Poisson impulse rate A = 1 /380T. We chose a large inter-impulse 
interval (380T) because the SL oscillators have a very slow relaxation back to the limit-cycle orbit following a perturbation. 



